package edu.amm.algorithms;

/**
 * Метод Холецкого решения СЛАУ.
 * 
 * @author Иванов Илья
 * @since 2013-05-12
 *
 */
public class Cholesky {
	
	private int dim;
	
	private double[][] A;
	private double[] b;
	
	private double[][] L;
	private double[][] U;
	
	/**
	 * Решает СЛАУ с заданной матрицей и вектором правой части
	 * методом Холецкого.
	 * 
	 * @param matrix Матрица системы.
	 * @param rightHand Вектор правой части системы.
	 * @return Решение системы.
	 */
	public double[] solve(double[][] matrix, double[] rightHand) {
		if (matrix == null || rightHand == null || matrix.length != rightHand.length)
			throw new RuntimeException("Некорректные параметры метода решения СЛАУ");
		
		checkMatrix(matrix);
		
		A = matrix;
		b = rightHand;
		dim = matrix.length;
		
		luDecomposition();
		double[] y = forwardTrace();
		double[] x = backwardTrace(y);
		
		return x;
	}
	
	private void checkMatrix(double[][] matrix) {
		int dim = matrix.length;
		for (int i = 0; i < dim; i++) {
			double[] row = matrix[i];
			if (row == null || row.length != dim)
				throw new RuntimeException("Матрица системы имеет некорректную конфигурацию");
		}
	}
	
	private void luDecomposition() {
		L = new double[dim][dim];
		U = new double[dim][dim];
		
		for (int k = 0; k < dim; k++) {
			for (int i = k; i < dim; i++) {
				double sum = 0;
				for (int s = 0; s < k; s++)
					sum += L[i][s] * U[s][k];
				
				L[i][k] = A[i][k] - sum;
			}
			
			for (int j = k + 1; j < dim; j++) {
				double sum = 0;
				for (int s = 0; s < k; s++)
					sum += L[k][s] * U[s][j];
				
				U[k][j] = (A[k][j] - sum) / L[k][k];
			}
			
			U[k][k] = 1;
			
			for (int s = k + 1; s < dim; s++) {
				L[k][s] = 0;
				U[s][k] = 0;
			}
		}
	}
	
	private double[] forwardTrace() {
		double[] y = new double[dim];
		for (int i = 0; i < dim; i++) {
			double s = 0;
			for (int j = 0; j < i; j++)
				s += L[i][j] * y[j];
			
			y[i] = (b[i] - s) / L[i][i];
		}
		
		return y;
	}
	
	private double[] backwardTrace(double[] y) {
		double[] x = new double[dim];
		for (int i = dim - 1; i >= 0; i--) {
			double S = 0;
			for (int j = i + 1; j < dim; j++)
				S += U[i][j] * x[j];
			
			x[i] = y[i] - S;
		}
		
		return x;
	}
}